function [F,G,Fz,Fzp,Fnp,Gz,Gzp,Gnp, ... 
    GFzz,GFzpz,GFnpnp,GFznp,GFzpnp,GFzpzp, ... 
    GFzpzpzp,GFzpzpnp,GFzpzpz,GFzpzz,GFzpznp,GFzpnpnp,GFzzz,GFzznp,GFznpnp,GFnpnpnp, ...
    ns,z_vec,zp_vec,np_vec] = model(varargin)



% temporary path name addition
warning off all
%addpath(path,'C:\kwre_so\lieps\adj\');
%addpath(path,'C:\kwre_so\lieps\adj\sym\');
warning on all

disp_eqns=0;
if disp_eqns
    de='';
else
    de=';';
end

% Define symbolic variables
disp(' ');
disp('Defining symbolic variables...');
tic
def_syms
disp(toc);

% Write equations
tic
disp(' ');
disp('Setting up model equations...');
eqn=0;






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Start of  : SETTING UP EQUATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% preliminaries

% utils
d1U = '(C^-sig)';
d1U_f = '(C_f^-sig)';
d2U = '(-chi*(N^eta))';

% production functions
F = '((K^al)*((A*N)^(1-al)))';
F1 = '(al*(K^(al-1))*((A*N)^(1-al)))';
F1_f = '(al*(K_f^(al-1))*((A_f*N_f)^(1-al)))';
F2 = '((1-al)*(K^al)*((A*N)^-al))';


%% Equations to be put in F(.)

% eq (1)
eqn=eqn+1;
val='log(xistar) - taustar';
eval(['f',num2str(eqn),' = ',val,de]);

% eq (2)
eqn=eqn+1;
val='phi1*taustar - Ftilde';
eval(['f',num2str(eqn),' = ',val,de]);

% eq (3)
eqn=eqn+1;
val='xistar_ss*fss*xistar - bigxi';
eval(['f',num2str(eqn),' = ',val,de]);

% eq (4)
eqn=eqn+1;
val='phi0*xistar - biglam';
eval(['f',num2str(eqn),' = ',val,de]);

% eq (22)
eqn=eqn+1;
val='Q*symG(itarg)-itarg-xistar';
eval(['f',num2str(eqn),' = ',val,de]);

% eq (23)
eqn=eqn+1;
val='Q*symdG(itarg)-sym(1)';
eval(['f',num2str(eqn),' = ',val,de]);

% eq (24)
eqn=eqn+1;
val = 'log(itarg)+log(K)+Ftilde-log(bigI)';
eval(['f',num2str(eqn),' = ',val,de]);


% eq (26a)
eqn=eqn+1;
val=[F,'-Y'];
eval(['f',num2str(eqn),' = ',val,de]);

% eq (26b)
eqn=eqn+1;
val='bigI+C+K*bigxi-Y';
eval(['f',num2str(eqn),' = ',val,de]);

% eq (27)
eqn=eqn+1;
val=['A*',F2,'+(',d2U,'/',d1U,')'];
eval(['f',num2str(eqn),' = ',val,de]);

% eq (28)
eqn=eqn+1;
val=['(b*',d1U_f,'/',d1U,')*(',F1_f,'+(1-delta)*Q_f+biglam_f)-Q'];
eval(['f',num2str(eqn),' = ',val,de]);



%%%% Auxiliary equations (not needed for solutions) %%%%

% definition for hazard function
eqn=eqn+1;
val='haz-exp(Ftilde)';
eval(['f',num2str(eqn),' = ',val,de]);

% definition for I/K (investment-capital ratio)
eqn=eqn+1;
val='IKrat-(bigI/K)';
eval(['f',num2str(eqn),' = ',val,de]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% Equations to be put in G(.)
% eq (25)
eqn=eqn+1;
val='(1-delta)*K+symG(itarg)*exp(Ftilde)*K-K_f';
eval(['f',num2str(eqn),' = ',val,de]);

% Productivity
eqn=eqn+1;
val='A_f - rhoa*A - avol*Ashock';
eval(['f',num2str(eqn),' = ',val,de]);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  End of  : SETTING UP EQUATIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%% define vector of contrals (y) and states (x)

x_vec = [K A];
y_vec = [Y N C bigI Q xistar itarg haz IKrat taustar Ftilde bigxi biglam];

xp_vec = [K_f A_f];
yp_vec = [Y_f N_f C_f bigI_f Q_f xistar_f itarg_f haz_f IKrat_f taustar_f Ftilde_f bigxi_f biglam_f];


z_vec = [y_vec,x_vec];
zp_vec = [yp_vec,xp_vec];
np_vec = [Ashock];

ns = length(x_vec);
nz = length(z_vec);



%% Create function f
F=[];
for ii=1:eqn-2
    eval(['F=[F; f',num2str(ii),'];']);
end

G=[];
for ii=eqn-1:eqn
    eval(['G=[G; f',num2str(ii),'];']);
end

nf = numel(F) ;
ng = numel(G) ;

disp(toc) ;
disp(' ') ;
disp(['Number of variables: ',num2str(length(z_vec))]) ; 
disp(['Number of states: ',num2str(ns)]) ;



%% construct analytical derivatives
[Gz,Gzp,Gnp,Fz,Fzp,Fnp, ... 
    GFzz,GFzpz,GFnpnp,GFznp,GFzpnp,GFzpzp, ... 
    GFzpzpzp,GFzpzpnp,GFzpzpz,GFzpzz,GFzpznp,GFzpnpnp,GFzzz,GFzznp,GFznpnp,GFnpnpnp] ...
     = adm(F,G,z_vec,zp_vec,np_vec,varargin{:}) ;


% for Miao_Wang trick (elas)

Fz_trick = Fz;
Fz_trick(2,10) = (-(phi1^2)+phi1+phi1*phi2)*taustar;
Fz_trick(3,6) = (fss+xistar_ss*dfss)*xistar;
Fz_trick(4,6) = (fss)*xistar;

[Fzz]=ad_trick(F,G,Fz_trick,z_vec,zp_vec,np_vec) ;
GFzz(ng+1:ng+nf,:) = Fzz ;

end % function


















